function  [lambda,V] = invpow(A,alpha,X,epsilon,max1,show)
% [lambda,V] = invpow(A,alpha,X,epsilon,max1,show)
% The inverse power method is used to find the dominant eigenpair.
% A is an n x n matrix, input.
% X is the starting vector, input.
% alpha is the given shift, input.
% epsilon is the tolerance, input.
% max1 is the maximum number of iterations, input.
% lambda is the dominant eigenvalue, output.
% V is the dominant eigenvector, output.
if nargin==5, show = 0; end
[n n] = size(A);
A = A - alpha*eye(n);               % Form the matrix A - aIpha I
[A,row] = LUfact(A);
lambda = 0;
cnt = 0;
err = 1;
done = 0;
iterating = 1;
state = iterating;
while ((cnt<=max1)&(state==iterating))
  Y = LUsolv(A,X,row);
  [m j] = max(abs(Y));
  c1 = Y(j);                        % Find "largest" element of Y.
  dc = abs(lambda - c1);
  Y = (1/c1)*Y;                     % Perform scalar multiplication
  dv = norm(X-Y);
  err = max(dc,dv);
  X = Y;                            % Update vector  X
  lambda = c1;                      % Update scalar  lambda 
  state = done;
  if (err>epsilon),                 % Check for convergence
    state = iterating;
  end
  cnt = cnt+1;
  lamb = alpha + 1/c1;
  if show==1,
    home; if cnt==1, clc; end; 
    disp(['Inverse Power Method iteration No. ',int2str(cnt)]),...
    disp(''),disp('c1 = '),disp(c1),...
    disp('lambda = '),disp(lamb),...
    disp('Eigenvector = '),disp(X)
  end
end
lambda = alpha + 1/c1;
V = X;
